מעבדה 1 - עיבוד תמונה

Rom Hirsch and Yarom Swissa

imports models & functions:

import cv2 as cv
import sys
import numpy as np
from matplotlib import pyplot as plt
from matplotlib.colors import NoNorm
from PIL import Image
import numpy.matlib
import pandas as pd
import math
from matplotlib import cm
from matplotlib.ticker import LinearLocator, FormatStrFormatter
from scipy import signal
from skimage import io

#%% functions
def tableDf(df):
    fig, ax = plt.subplots()
    # hide axes
    fig.patch.set_visible(False)
    ax.axis('off')
    ax.axis('tight')
    ax.table(cellText=df.values, colLabels=df.columns, loc='center')
    fig.tight_layout()

def imageShow(im,  title = "image"):
    """display image"""
    cv.imshow(title, im)
    cv.waitKey(0)
    cv.destroyAllWindows()

def uint8c(x):
    """custom uint8 - (if val < 0 then val = 0)"""
    x2 = np.copy(x)
    x2[x2 < 0] = 0
    x2[x2 > 255] = 255
    return np.uint8(x2)



def plotPlusImage(axP, axI, axN, im, title):
    axP.plot(im[1, :])
    axP.set_title(title)
    axI.imshow(im, cmap="gray", norm=NoNorm())
    axI.set_title(title)
    imgMin = np.min(im.flatten())
    imgMax = np.max(im.flatten())
    normalImg = np.dot(255, (im - imgMin) / (imgMax - imgMin))
    axN.imshow(normalImg, cmap="gray")
    axN.set_title(f"normalize {title}")
    #axN.plot(normalImg[0, :])


# plot stem with color
def stem_plot(n, val, color):
    markerline1, stemlines1, baseline1 = plt.stem(n, val)
    plt.setp(markerline1, 'markerfacecolor', color)
    plt.setp(stemlines1, linestyle="-", color=color, linewidth=2)

# Create Stem plot with color
def plotStem(title, ylabel, xlabel, color, x, y):
    plt.title(title)
    plt.ylabel(ylabel)
    plt.xlabel(xlabel)
    stem_plot(x, y, color)
    #plt.show()

#padded zeros to array , params (array ,size of the array after the padded)
def padded_zeros(arr, N): # padded zeros
    padded_array = np.zeros(N)
    padded_array[:len(arr)] += arr
    return padded_array

def add_sp(img, probability = 0.1):
  noise = np.random.uniform(size = img.shape)
  salt = noise > 1 - 0.5 * probability
  pepper = noise < 0.5 * probability
  return np.clip(img.astype('float64') + salt * 255 - pepper * 255, 0, 255).astype('uint8')
  salt_pepper_noisy = add_sp(image)
  
def random_gaussian(N = 100000, mu = 0, sigma = 1):
  x = np.zeros(N)
  for i in range(0, N):
    x[i] += (np.random.normal(0, 1) * sigma) + mu
  return x

def conv2d(image, filter, mode = 'full', add_offset = 0, clip_enabled = False, clip_min = 0, clip_max = 255):
  pad_amnt_y = len(filter) - 1
  pad_amnt_x = len(filter[0]) - 1
  padded_img = np.zeros((len(image) + 2 * pad_amnt_y, len(image[0]) + 2 * pad_amnt_x))
  padded_img[pad_amnt_y : len(padded_img) - pad_amnt_y, pad_amnt_x : len(padded_img[0]) - pad_amnt_x] = image.astype('float64')
  output = np.zeros((len(image) + pad_amnt_y * 2, len(image[0]) + pad_amnt_x * 2))
  for x in range(0, len(output[0]) - len(filter[0]) + 1):
    for y in range(0, len(output) - len(filter) + 1):
      output[y][x] = np.sum(padded_img[y: y + len(filter), x : x + len(filter[0])] * filter)
  output_float = output[:len(output) - pad_amnt_y, :len(output[0]) - pad_amnt_x]
  output = output_float + add_offset
  if (clip_enabled):
    output = np.clip(output, clip_min, clip_max)
  return output.astype('uint8'), output_float

1 משימת חימום

1.1 בניית תמונה - ריבוע בתוך ריבוע

def SquareInSquare(color_backgroud, color_object, res = False, range_background = 400):
    black_with_square = np.ones([range_background, range_background]) * color_backgroud
    black_with_square[150:250, 150:250] = color_object
    black_with_square = np.uint8(black_with_square)
    if res:
        Cw = (color_object - color_backgroud) / color_backgroud
        return black_with_square, Cw
    else:
        fig = plt.figure(figsize=(10, 10))
        plt.subplot(211)
        plt.imshow(black_with_square, cmap="gray", norm=NoNorm())
        plt.title(f"SquareInSquare Bg = {color_backgroud} Bo = {color_object}")
        plt.subplot(212)
        vals = black_with_square.flatten()
        plt.hist(vals, 256, [0, 256])
        Cw = (color_object - color_backgroud) / color_backgroud
        plt.title(f"Waber Contrast is:{Cw}")
        #plt.show()

SquareInSquare(220,50)

1.2 ייצור תמונת פסים

num_stripes = np.arange(0, 256, int(256/8))
        print(num_stripes)

[  0  32  64  96 128 160 192 224]

def stripes_uint8(height, width, num_stripes):
    if np.isscalar(num_stripes):
        num_stripes = np.arange(0, 256, int(256/num_stripes))
    width_stripe = math.floor(width/len(num_stripes))
    vector_row = np.ones(width)
    vector_row = np.array([vector_row[(i * width_stripe): ((i + 1) * width_stripe)] * v for i, v in enumerate(num_stripes)])
    vector_row = vector_row.flatten()
    return np.uint8(np.matlib.repmat(vector_row,height,1))


#1.
height = 200
width = 300
num_stripes = [0, 50, 100, 150, 200, 250]
plt.figure()
plt.imshow(stripes_uint8(height, width, num_stripes), cmap="gray")
plt.title("Stripes")

Text(0.5, 1.0, 'Stripes')

1.3 .יצירה והצגה של אותות

1.3.1 הצגת אות סינוס דו מימד בפורמטים וטיפוסי משתנים שונים

N = 256
Fs = N
Ts = 1/Fs
fx = 5
t = np.linspace(0, N-1, N) * Ts
sx = N
imX = np.zeros((N, N))
imX = np.matlib.repmat(50*np.sin(2*np.pi*fx*t), N, 1)

# plot
fig, axs = plt.subplots(3, 4,figsize=(16,13))
plotPlusImage(axs[0, 0], axs[1, 0], axs[2, 0], uint8c(np.abs(imX)), 'A line of uint(abs(imX))')
plotPlusImage(axs[0, 1], axs[1, 1], axs[2, 1], uint8c(imX), 'A line of uint(imX)')
plotPlusImage(axs[0, 2], axs[1, 2], axs[2, 2], np.float32(imX), 'A line of imX')
plotPlusImage(axs[0, 3], axs[1, 3], axs[2, 3], uint8c(imX-np.min(imX.flatten())), 'A line of uint(imX-min(imX.flatten()))')
plt.show()

1.3.2 יצירה והצגה של אותות חד מימדיים בזמן ובתדר

# question 1 -4 
A = 7
f = 10
L = 1
fs = 200
Ts = 1/fs
t = np.arange(0,L-Ts,Ts)
#t = np.linspace(0, L, int(L/Ts))
y = np.cos(2*np.pi*f*t)
N = len(t)
plt.figure(figsize=(15, 10))
plt.subplot(221)
plt.plot(t,y)
plt.title(f"y = sin(2*pi*f*t) fs = {fs}")
plt.subplot(222)
plt.plot(np.linspace(-(fs/2), fs/2-1, y.shape[0]), np.fft.fftshift(abs(np.fft.fft(y))))
plt.title(f"(abs(fft(y)) fs = {fs}")

fs = 12
Ts = 1/fs
t = np.arange(0,L-Ts,Ts)
y = np.cos(2*np.pi*f*t)
plt.subplot(223)
plt.plot(t, y)
plt.title(f"y = sin(2*pi*f*t) fs = {fs}")
plt.ylim(-1,1)
plt.subplot(224)
plt.plot(np.linspace(-(fs/2), fs/2-1, y.shape[0]), np.fft.fftshift(abs(np.fft.fft(y))))
plt.title(f"(abs(fft(y)) fs = {fs}")

Text(0.5, 1.0, '(abs(fft(y)) fs = 12')

אות בטווח של 0-255 + unit8

#%% 5
A = 7
f = 10
L = 1
fs = 200
Ts = 1/fs
t = np.arange(0,L-Ts,Ts)
y = np.cos(2*np.pi*f*t)
y = uint8c((1 + y) * 128)
N = len(t)
plt.figure(figsize=(10, 10))
plt.subplot(221)
plt.plot(t,y)
plt.title(f"y = unit8(sin(2*pi*f*t)) fs = {fs}")
plt.subplot(222)
plt.plot(np.linspace(-(fs/2), fs/2-1, y.shape[0]), np.fft.fftshift(abs(np.fft.fft(y))))
plt.title(f"(abs(fft(y)) fs = {fs}")
fs = 12
Ts = 1/fs
t = np.arange(0,L-Ts,Ts)
y = np.cos(2*np.pi*f*t)
y = uint8c((1 + y) * 128)
N = len(t)
plt.subplot(223)
plt.plot(t,y)
plt.title(f"y = unit8(sin(2*pi*f*t)) fs = {fs}")
plt.subplot(224)
plt.plot(np.linspace(-(fs/2), fs/2-1, y.shape[0]), np.fft.fftshift(abs(np.fft.fft(y))))
plt.title(f"(abs(fft(y)) fs = {fs}")

Text(0.5, 1.0, '(abs(fft(y)) fs = 12')

6

רק האות בתדירות דגימה של 200 הרץ עומד בתנאי נייקוויסט שהוא דורש תדר כפול מתדר ההאות כלומר המינימום שלא יווצרו לנו התחזויות הוא אות תגימה של 20

בדגימה של 200 קיבלנו בספקטרות תדר 10 כפי שאות המקור לאומת זאת בדגימה של 12 קיבלנו אות בתדר 2 זה נובע בגלל שהפרש בין התדגימה לאות הוא 2 לכן ציפינו לקבל אות בתדירות 2

בשאלות 5 כפי שציפינו נוסף רכיב דיסי לספקטרום בגלל הנירמול בין 0-255

1.4 קונבולוציה חד מימדית, אות, רעש וסיכון חד מימדיים

1.4.1 - נוסחת הקונבולוציה

1.4.2 - קונבולוציה חד מימדית

2 - פיתרון ידני של הקונבולוציה

3 - פיתרון ע"י קונבולציה גרפית

שרטוט גרפית של האותות והתוצאה:

f = np.array([1, 2, 3, 1, 2])
h = np.array([1, 1])
g = np.convolve(f, h)
n = np.arange(0, 4, 1)
m1 = len(f)
m2 = len(h)
f = padded_zeros(f, m2+m1-1)
h = padded_zeros(h, m2+m1-1)
plt.figure(figsize=(10, 10))
plt.subplot(311)
plotStem("f[n]", "amplitude", "n", 'red', range(m1+m2-1), f)
plt.subplot(312)
plotStem("h[n]", "amplitude", "n", 'green', range(m1+m2-1), h)
convXY = np.convolve(f, h)
plt.subplot(313)
plotStem("conv(f*h)", "amplitude", "n", 'blue', range(convXY.shape[0]), convXY)
plt.show()

/usr/local/lib/python3.7/dist-packages/ipykernel_launcher.py:41: UserWarning: In Matplotlib 3.3 individual lines on a stem plot will be added as a LineCollection instead of individual lines. This significantly improves the performance of a stem plot. To remove this warning and switch to the new behaviour, set the "use_line_collection" keyword argument to True.
/usr/local/lib/python3.7/dist-packages/ipykernel_launcher.py:41: UserWarning: In Matplotlib 3.3 individual lines on a stem plot will be added as a LineCollection instead of individual lines. This significantly improves the performance of a stem plot. To remove this warning and switch to the new behaviour, set the "use_line_collection" keyword argument to True.
/usr/local/lib/python3.7/dist-packages/ipykernel_launcher.py:41: UserWarning: In Matplotlib 3.3 individual lines on a stem plot will be added as a LineCollection instead of individual lines. This significantly improves the performance of a stem plot. To remove this warning and switch to the new behaviour, set the "use_line_collection" keyword argument to True.

5-6 קונבולוציה עם הלם מוזז

# 1.4 - q6
f = np.array([1, 2, 3, 1, 2])
h2 = np.zeros(f.shape[0]*2)
m1 = len(f)
m2 = len(h2)
h2[[0,f.shape[0]]] = 1
f = padded_zeros(f, m2)
plt.figure(figsize=(10, 10))
plt.subplot(311)
plotStem("f[n]", "amplitude", "n", 'red', range(m2), f)
plt.subplot(312)
plotStem("h[n]", "amplitude", "n", 'green', range(m2), h2)
g2 = np.convolve(f,h2)
plt.subplot(313)
plotStem("conv(f*h)", "amplitude", "n", 'blue', range(m2), g2[0:m2])
plt.show()

/usr/local/lib/python3.7/dist-packages/ipykernel_launcher.py:41: UserWarning: In Matplotlib 3.3 individual lines on a stem plot will be added as a LineCollection instead of individual lines. This significantly improves the performance of a stem plot. To remove this warning and switch to the new behaviour, set the "use_line_collection" keyword argument to True.
/usr/local/lib/python3.7/dist-packages/ipykernel_launcher.py:41: UserWarning: In Matplotlib 3.3 individual lines on a stem plot will be added as a LineCollection instead of individual lines. This significantly improves the performance of a stem plot. To remove this warning and switch to the new behaviour, set the "use_line_collection" keyword argument to True.
/usr/local/lib/python3.7/dist-packages/ipykernel_launcher.py:41: UserWarning: In Matplotlib 3.3 individual lines on a stem plot will be added as a LineCollection instead of individual lines. This significantly improves the performance of a stem plot. To remove this warning and switch to the new behaviour, set the "use_line_collection" keyword argument to True.

1.4.3 יצירת אות חד ממדי (רעש) בעל התפלגות נתונה (גאוסית לבנה)

# Gaussian Noise
mu = np.array([0, 3, 5, 10])
sigma = np.array([0.5, 2.5, 9, 10])
noise = np.array([np.random.normal(mu[i], sigma[i], 1000) for i in range(mu.shape[0])])
plt.figure(figsize=(14, 14))
for i in range(noise.shape[0]):
    plt.subplot(4, 1, i+1)
    plt.plot(noise[i])
    plt.title(f"noise {i}")
    plt.xlabel("samples", horizontalalignment='right', x=1.0)
    plt.ylabel("amp")
plt.show()

plt.figure(figsize=(15, 15))
for i in range(noise.shape[0]):
    plt.subplot(4, 1, i+1)
    counts, bins = np.histogram(noise[i])
    plt.hist(bins[:-1], bins, weights=counts)
    plt.title(f"noise {i}")
    plt.xlabel("samples", horizontalalignment='right', x=1.0)
    plt.ylabel("amp")
plt.show()

1.4.4 סינון אות חד ממדי רועש

8 הוספת רעש גאוסי והורדת הרעש ע"י מסנן

#%% 1.4.4

h1 = np.array([0.2]*5)
A = 7
f = 10
L = 1
fs = 200
Ts = 1/fs
t = np.arange(0,L-Ts,Ts)
y = np.cos(2*np.pi*f*t)
n1 = np.random.normal(mu[0], sigma[0], len(y))
f1 = n1 + y 
g1 = np.convolve(f1,h1)
plt.figure(figsize=(20,15))
plt.subplot(413)
plt.plot(f1)
plt.title("f1[n]")
plt.subplot(411)
plt.plot(n1)
plt.title("n1[n]")
plt.subplot(412)
plt.plot(y)
plt.title("y[n]")
plt.subplot(414)
plt.plot(g1)
plt.title("g1[n]")

Text(0.5, 1.0, 'g1[n]')

9 - מסקנות

מסקנות : הרעש שהתווסף הוא רעש גאוסי אשר מתפרס על כל התדרים ולכן לאחד שהעברנו את האות במסנן מעביר נמוכים נפתרנו מרעש בתדרים הגבוהים אך עדיין נשאר הרעש בתדרים הנמוכים

ממדי האות לאחר הקונבולוציה יצאו הגודל של אות המקורי ועוד הגודל של המסנן פחות אחד וזה מתאים לתאוריה

plt.figure(figsize=(20,15))
plt.subplot(311)
plt.title("h1[n]")
plt.plot(np.linspace(-fs/2, fs/2,h1.shape[0]),np.fft.fftshift(np.fft.fft(h1)))
plt.subplot(312)
plt.title("f1[n]")
plt.plot(np.linspace(-fs/2, fs/2,f1.shape[0]),np.fft.fftshift(np.fft.fft(f1)))
plt.subplot(313)
plt.title("g1[n]")
plt.plot(np.linspace(-fs/2, fs/2,g1.shape[0]),np.fft.fftshift(np.fft.fft(g1)))

/usr/local/lib/python3.7/dist-packages/matplotlib/cbook/__init__.py:1317: ComplexWarning: Casting complex values to real discards the imaginary part
  return np.asarray(x, float)
/usr/local/lib/python3.7/dist-packages/matplotlib/cbook/__init__.py:1317: ComplexWarning: Casting complex values to real discards the imaginary part
  return np.asarray(x, float)
/usr/local/lib/python3.7/dist-packages/matplotlib/cbook/__init__.py:1317: ComplexWarning: Casting complex values to real discards the imaginary part
  return np.asarray(x, float)

[<matplotlib.lines.Line2D at 0x7f7a4ddd2ad0>]

2 הניסוי

2.1 בהירות,ניגודיות ומערכת הראייה האנושית

קוד -המששמש לאפליקציה אינטרקטיבית לקביעת ניגודיות וובר

interactive application :

control keys:

object brightness: z - increase, x - decrease ;

background brightness: a - increase, s - decrease, g - display histogram + save state to the table q - exit

#%% interactive application : key: object z - increase, x - decrease ; 
# background  a - increase, s - decrease, 
# g - display histogram + save state to the table, q - exit

bg=1
bo=1  
d = {"brightness (bg)" : ["brightness (bo)", "JND(Cw)"]}
c = 1
# to 
"""
while (True):
    image,Cw = SquareInSquare(bg,bo,True)
    cv.imshow(f"bg = {bg}, bo = {bo}",image)   
    k = cv.waitKey(0)
    if k == ord('a'): 
        bg += 1
    elif k == ord('s'): 
        bg = bg - 1
    elif k == ord('z'): 
        bo += 1
    elif k == ord('x'): 
        bo = bo - 1
    elif k == ord('q'):
        cv.destroyAllWindows()
        break
    elif k == ord('g'):
        #addToTable()
        if f"{bg}" in d.keys(): 
            d[f"{bg}({c})"] = [bo,Cw]
        else:
            d[f"{bg}"] = [bo,Cw]
        SquareInSquare(bg,bo)
    else:
        continue
    cv.destroyAllWindows()
df = pd.DataFrame(data=d)
print(df.to_markdown()) # require pip install tabulate 
"""
SquareInSquare(220,180)

התוצאות שלנו המניסוי: במעבדה בכיתה באור דלוק

#%%resulte of the test: 
d = {"brightness (bg)" : ["brightness (bo)", "JND(Cw)"], '40': [39, -0.025],'40(1)': [41, 0.025],'140': [140, 0.0],'140(1)': [141, 0.007142857142857143],'220': [221, 0.004545454545454545],'220(1)': [220, 0.0]}
df = pd.DataFrame(data=d, index=['brightness(bo)', 'Cw'])
print(df.to_markdown(index=False)) # require pip install tabulate

| brightness (bg)   |     40 |   40(1) |   140 |       140(1) |          220 |   220(1) |
|:------------------|-------:|--------:|------:|-------------:|-------------:|---------:|
| brightness (bo)   | 39     |  41     |   140 | 141          | 221          |      220 |
| JND(Cw)           | -0.025 |   0.025 |     0 |   0.00714286 |   0.00454545 |        0 |

תוצאות מבדיקה נוספת בבית בתנאים אחרים:

d = {'brightness (bg)': ['brightness (bo)', 'JND(Cw)'],
 '40': [41, 0.025],
 '140': [139, -0.007142857142857143],
 '140(1)': [142, 0.014285714285714285],
 '220': [219, -0.004545454545454545],
 '220(1)': [221, 0.004545454545454545]}
df = pd.DataFrame(data=d, index=['brightness(bo)', 'Cw'])
print(df.to_markdown(index=False)) # require pip install tabulate

| brightness (bg)   |     40 |          140 |      140(1) |          220 |       220(1) |
|:------------------|-------:|-------------:|------------:|-------------:|-------------:|
| brightness (bo)   | 41     | 139          | 142         | 219          | 221          |
| JND(Cw)           |  0.025 |  -0.00714286 |   0.0142857 |  -0.00454545 |   0.00454545 |

גרף תאורתי של רגישות :

התדר הדומיננטי בגל ריבועי הוא תדר מאוד נמוך לכן הC

The dominant frequency of square wave is near DC, so by looking at the theoretical graph we would expect CW around 30 A

כושר ההבחנה תלוי גם בתאורה בחדר בסוג המסך ובצופה לכן יתקבל תוצאות שונות לכל בדיקה במידה ותנאים אלו ישתנו

2.1.2 המחשת תכונת הגזירה של העין האנושית והבחנה בין רמות אפות פסי מאך

f, ax = plt.subplots(4, 3,figsize=(25,25))

img=stripes_uint8(height, width, 8)
ax[0, 0].imshow(img,cmap="gray")
ax[0, 0].set_title("8 bands Image")
ax[0, 0].axis('off')
print(img.flatten()[-1])
ax[0, 1].hist(img.flatten(), 256, [0, 256])
ax[0, 1].set_title("8 bands Histogram")
ax[0, 2].plot(img[0,:])
ax[0, 2].set_title("the 1st line")

img=stripes_uint8(height,width,16)
ax[1, 0].imshow(img,cmap="gray")
ax[1, 0].set_title("16 bands Image")
ax[1, 0].axis('off')
ax[1, 1].hist(img.flatten(), 256, [0, 256])
ax[1, 1].set_title("8 bands Histogram")
ax[1, 2].plot(img[0,:])
ax[1, 2].set_title("the 1st line")

img=stripes_uint8(height, width, 32)
ax[2, 0].imshow(img,cmap="gray")
ax[2, 0].set_title("32 bands Image")
ax[2, 0].axis('off')
ax[2, 1].hist(img.flatten(), 256, [0, 256])
ax[2, 1].set_title("32 bands Histogram")
ax[2, 2].plot(img[0,:])
ax[2, 2].set_title("the 1st line")

img=stripes_uint8(height, width, 256)
ax[3, 0].imshow(img,cmap="gray")
ax[3, 0].set_title("Ramp Image")
ax[3, 0].axis('off')
ax[3, 1].hist(img.flatten(), 256, [0, 256])
ax[3, 1].set_title("Ramp Histogram")
ax[3, 2].plot(img[0,:])
ax[3, 2].set_title("the 1st line")
plt.show()
img=stripes_uint8(height, width, 8) 
noise = np.uint8(np.matlib.repmat(np.random.normal(0, 20, img.shape[0]),img.shape[1],1))

224

הוספת רעש:

img=stripes_uint8(height, width, 8)
row,col= img.shape
mean = 0
sigma = 20
gauss = np.random.normal(mean,sigma,(row,col))
gauss = gauss.reshape(row,col)



f, ax = plt.subplots(4, 3,figsize=(25,25))

img=stripes_uint8(height, width, 8)
img = img + gauss
ax[0, 0].imshow(img,cmap="gray")
ax[0, 0].set_title("8 bands Image")
ax[0, 0].axis('off')
print(img.flatten()[-1])
ax[0, 1].hist(img.flatten(), 256, [0, 256])
ax[0, 1].set_title("8 bands Histogram")
ax[0, 2].plot(img[0,:] )
ax[0, 2].set_title("the 1st line")

img=stripes_uint8(height,width,16)
row,col= img.shape
mean = 0
sigma = 20
gauss = np.random.normal(mean,sigma,(row,col))
gauss = gauss.reshape(row,col)
img = img + gauss
ax[1, 0].imshow(img,cmap="gray")
ax[1, 0].set_title("16 bands Image")
ax[1, 0].axis('off')
ax[1, 1].hist(img.flatten(), 256, [0, 256])
ax[1, 1].set_title("8 bands Histogram")
ax[1, 2].plot(img[0,:])
ax[1, 2].set_title("the 1st line")

img=stripes_uint8(height, width, 32)
row,col= img.shape
mean = 0
sigma = 20
gauss = np.random.normal(mean,sigma,(row,col))
gauss = gauss.reshape(row,col)
img = img + gauss
ax[2, 0].imshow(img,cmap="gray")
ax[2, 0].set_title("32 bands Image")
ax[2, 0].axis('off')
ax[2, 1].hist(img.flatten(), 256, [0, 256])
ax[2, 1].set_title("32 bands Histogram")
ax[2, 2].plot(img[0,:])
ax[2, 2].set_title("the 1st line")

img=stripes_uint8(height, width, 256)
row,col= img.shape
mean = 0
sigma = 20
gauss = np.random.normal(mean,sigma,(row,col))
gauss = gauss.reshape(row,col)
img = img + gauss
ax[3, 0].imshow(img,cmap="gray")
ax[3, 0].set_title("Ramp Image")
ax[3, 0].axis('off')
ax[3, 1].hist(img.flatten(), 256, [0, 256])
ax[3, 1].set_title("Ramp Histogram")
ax[3, 2].plot(img[0,:])
ax[3, 2].set_title("the 1st line")
plt.show()
img=stripes_uint8(height, width, 8) 
noise = np.uint8(np.matlib.repmat(np.random.normal(0, 20, img.shape[0]),img.shape[1],1))

266.1002506128287

2.2 אותות דו מימדיים, ומערכת הראייה

2.2.1 יצירת הצגה של אותות דו

class d2signal:
    def __init__(self, sx, sy, fs, fx, fy, amp=1, phase=0, func=np.sin):
        self.n = N
        self.amp = amp
        self.fx = fx
        self.fy = fy
        self.fs = fs
        self.func = func
        self.sx = np.array(np.linspace(0, sx-1, sx)).reshape((1,-1))
        self.sy = np.array(np.linspace(0, sy-1, sy).reshape((-1,1))) 
        self.img = self._genSignal()
        
    def _genSignal(self):
        Ts = 1/self.fs
        wx = 2 * np.pi * self.fx * Ts;
        wy = 2 * np.pi * self.fy * Ts;        
        self.img =  127 -  self.amp * self.func(wy * self.sy + wx * self.sx)
        self._signalUint8()
        return self.img
    
    def _signalUint8(self):
        self.imgUint8 = uint8c(self.img) 

    def imshow(self,title="image", Uint8 = False):
        if Uint8:
            img = self.imgUint8
        else:
            img = self.img
        plt.figure(figsize=(15,15))
        plt.title(title)
        plt.imshow(img, cmap="gray",norm=NoNorm())
        plt.show()
        
    def normalize(self):
        imgMin = np.min(self.img.flatten())
        imgMax = np.max(self.img.flatten())
        self.old = self.img.copy()
        self.img = np.dot(255, (self.img - imgMin) / (imgMax - imgMin))
        self._signalUint8()
        
        
    def fft(self, title="fft", Uint8 = False):
        plt.figure(figsize=(15,15))
        if Uint8:
            img = self.imgUint8
        else:
            img = self.img
        plt.title(title)
        plt.imshow(np.fft.fftshift(abs(np.fft.fft2(img))), cmap="gray")
        plt.show()
        
    def imshow3d(self, Uint8 = False):
        from matplotlib import cm
        from matplotlib.ticker import LinearLocator
        if Uint8:
            img = self.imgUint8
        else:
            img = self.img
        xx, yy = np.mgrid[0:img.shape[0], 0:img.shape[1]] # get an x y thing of the image
        
        fig = plt.figure(figsize=(15,15))
        ax = fig.gca(projection='3d') 
        surf = ax.plot_surface(xx, yy, img, cmap=cm.gray,
                               linewidth=0, antialiased=False)
        ax.set_zlim(0, 255)
        ax.zaxis.set_major_locator(LinearLocator(10)) 
        ax.zaxis.set_major_formatter(FormatStrFormatter('%d')) 
        fig.tight_layout() 
        plt.show()


I1 = d2signal(256,256,256,5,0,50)
I1.imshow("Amp = 50 , fx = 5",True)

#%%
I2 = d2signal(256,256,256,0,10,100)
I2.imshow("Amp = 100 fy= 10",True)
#%%
I3 = d2signal(256,256,256,5,10,127)
I3.imshow("Amp = 127 fy= 10 fx = 5",True)

מה שניתן לראות בתמונה 3 זה שיש לנו תדר שהוא בשני הצירים. לפי הספקטרום התדר ניתן לקבל את התדר המרחבי ואת הרכיבים בכל ציר

I1.fft()
I2.fft()
I3.fft()

I1.normalize()  
I2.normalize()  
I3.normalize() 
I1.imshow("Amp = 50 , fx = 5",True)
I2.imshow("Amp = 100 fy= 10",True)
I3.imshow("Amp = 127 fy= 10 fx = 5",True)

I4 = d2signal(256,256,256,5,0,10)
I4.img = I1.img + I2.img + I3.img  
I4.normalize()
I4.imshow("Amp = 50 , fx = 5",True)
I4.fft()

I4.imshow3d()

2.2.2 Campbell Robson

min_freq = 0.1 # cycles per image
N=400
min_amplitude = 1
max_amplitude = 127
max_freq = 40 
f = np.exp(np.linspace(np.log(min_freq), np.log(max_freq), N).reshape((1,-1)))
a = np.exp(np.linspace(np.log(min_amplitude), np.log(max_amplitude), N).reshape((-1,1)))
img = 128 + a * np.sin(2 * np.pi * f)
plt.figure(figsize=(15,15))
plt.imshow(img, cmap='gray', vmin=0, vmax=255)

<matplotlib.image.AxesImage at 0x7f7a4fbb1290>

2.3 2D Convolution

2.3.3 - Filters

filters + mse + snr + psnr:

url = r"https://placekitten.com/500/300"
image = io.imread(url,0) 
image = cv.cvtColor(image, cv.COLOR_BGR2GRAY)
#create filters
lpf = np.array([[1,1,1],[1,1,1],[1,1,1]]) / 9.0
hpf = np.array([[-1,-1,-1],[-1,8,-1],[-1,-1,-1]]) / 9.0
sharp_f = np.array([[-1,-1,-1],[-1,17,-1],[-1,-1,-1]]) / 9.0
laplacian_f = np.array([[0,-1,0],[-1,4,-1],[0,-1,0]])
im_lp, _ = conv2d(image, lpf)
im_laplacian, im_laplacian_float = conv2d(image, laplacian_f, add_offset = 127, clip_enabled = True, clip_min = 0, clip_max = 255)
im_hp, im_hp_float = conv2d(image, hpf, add_offset = 127, clip_enabled = True, clip_min = 0, clip_max = 255)
median_filt = cv.medianBlur(image, 3)
salt_pepper_noisy = add_sp(image, 0.1)
gaussian_noise = random_gaussian(N = (len(image) * len(image[0])), sigma = 25)
gaussian_noise = np.reshape(gaussian_noise, image.shape)
gaussian_noisy = np.clip((image.astype('float64') + gaussian_noise), 0, 255).astype('uint8')
salt_pepper_lpf, _ = conv2d(salt_pepper_noisy, lpf)
gaussian_lpf, _ = conv2d(gaussian_noisy, lpf)
salt_pepper_hpf, salt_pepper_hpf_float = conv2d(salt_pepper_noisy, hpf, add_offset = 127, clip_enabled = True, clip_min = 0, clip_max = 255)
gaussian_hpf, gaussian_hpf_float = conv2d(gaussian_noisy, hpf, add_offset = 127, clip_enabled = True, clip_min = 0, clip_max = 255)
gaussian_sharp, _ = conv2d(gaussian_noisy, sharp_f, clip_enabled = True)
salt_pepper_sharp, _ = conv2d(salt_pepper_noisy, sharp_f, clip_enabled = True)
salt_pepper_laplacian, salt_pepper_laplacian_float = conv2d(salt_pepper_noisy, laplacian_f, add_offset = 127, clip_enabled = True, clip_min = 0, clip_max = 255)
gaussian_laplacian, gaussian_laplacian_float = conv2d(gaussian_noisy, laplacian_f, add_offset = 127, clip_enabled = True, clip_min = 0, clip_max = 255)
salt_pepper_median = cv.medianBlur(salt_pepper_noisy, 3)
im_sharp, _ = conv2d(image, sharp_f, clip_enabled = True, clip_min = 0, clip_max = 255)

gaussian_median = cv.medianBlur(gaussian_noisy, 3)
salt_pepper_fixed, _ = conv2d(salt_pepper_median, sharp_f, clip_enabled = True)

# calc mse + snr + snr
img_ref = image
reference_max = np.max(img_ref).astype('float64') # 
reference_avg = np.average(img_ref).astype('float64')
reference_variance = np.sum((image.astype('float64') - reference_avg) * (image.astype('float64') - reference_avg) / (len(image) * len(image[0])))
plt.figure(figsize=(15,15))
figure, axes = plt.subplots(3, 3,figsize=(15,15))
axes[0][0].imshow(img_ref, cmap = 'gray', vmin = 0, vmax = 255)
axes[0][0].set(title = 'original')
axes_index = 1
imgs = [im_lp, im_hp_float, im_sharp, im_laplacian_float, gaussian_noisy, gaussian_lpf, salt_pepper_noisy, salt_pepper_fixed]
labels = ['orig - lpf', 'orig - hpf', 'orig - sharp', 'orig - laplacian', 'gaussian noise', 'gaussian - lpf', 'sp - noise', 'sp - nmedian&nsharp']
snr_log = []
psnr_log = []
mse = []
for img in imgs:
  # because the convolution increased the size of the image, ignore the borders
  if ((len(img) != len(img_ref)) and (len(img[0]) != len(img_ref[0]))):
    img_ = img[1:len(img) - 1,1:len(img[0]) - 1]
  else:
    img_ = img
  noise = img_ref.astype('float64') - img_.astype('float64')
  noise_avg = np.average(noise)
  noise_variance = np.sum((noise - noise_avg) * (noise - noise_avg)) / (len(noise) * len(noise[0]))
  snr_log.append(10*np.log10(reference_variance / noise_variance))
  current_mse = np.sum(noise * noise) / (len(noise) * len(noise[0]))
  mse.append(current_mse)
  psnr_log.append(10*np.log10(reference_max * reference_max / current_mse))
  axes[int(axes_index / 3)][int(axes_index % 3)].imshow(img_, cmap = 'gray', vmin = 0, vmax = 255)
  axes[int(axes_index / 3)][int(axes_index % 3)].set(title = labels[axes_index - 1])
  axes_index += 1
figure.tight_layout()

<Figure size 1080x1080 with 0 Axes>

כאשר מופיע רעש גאוסייני סדר הפעולות הסינון היעיל ביותר יהיה מסנן מעביר נמוכים ואז מסנן sharp

כאשר מופיע רעש מלח פלפל היעיל ביותר יהיה פילטר חציון

plt.figure(figsize=(15,10))
plt.xticks(np.arange(len(labels)),labels)
plt.plot(snr_log)
plt.plot(psnr_log)
plt.title('SNR, PSNR, log10 scale')
plt.legend(['SNR','PSNR'])
plt.show()

plt.figure(figsize=(15,10))
plt.xticks(np.arange(len(labels)),labels)
plt.plot(mse)
plt.title("MSE")
plt.show()

3 דו"ח מכין לניסוי והשלמות קוד

3.2

מספר הרמות המצג בתרשים פסי מאך ללא רעש מוגדר כך שככל שמספר רמות האפור גדול יותר העין מתפקדת כאינטגרטור ולכן יהיה יותר קשה להבחין בין פס לפס

לעומת זאת ככל שמספר רמות האפור קטן יותר העין מתפקדת הגוזר וניתן לראות את המעברים בין פס לפס בקלות יותר. העין יכולה להבחין בערך 100 רמות אפור

3.3

b.

A = 1
f = 66.67
L = 1
fs = 200
Ts = 1/fs
t = np.arange(0,L-Ts,Ts)
#t = np.linspace(0, L, int(L/Ts))
y = A * np.cos(2*np.pi*f*t)
N = len(t)
print(f"min: {np.min(y)} max: {np.max(y)}")
plt.figure(figsize=(15, 10))
plt.subplot(221)
plt.plot(t,y)
plt.title(f"y = sin(2*pi*f*t) fs = {fs}")
plt.subplot(222)
plt.plot(np.linspace(-(fs/2), fs/2-1, y.shape[0]), np.fft.fftshift(abs(np.fft.fft(y))))
plt.title(f"(abs(fft(y)) fs = {fs}")

min: -0.517668669793734 max: 1.0

Text(0.5, 1.0, '(abs(fft(y)) fs = 200')

c.

A = 1
f = 65.93
L = 1
fs = 200
Ts = 1/fs
t = np.arange(0,L-Ts,Ts)
#t = np.linspace(0, L, int(L/Ts))
y = A * np.cos(2*np.pi*f*t)
N = len(t)
print(f"min: {np.min(y)} max: {np.max(y)}")
plt.figure(figsize=(15, 10))
plt.subplot(221)
plt.plot(t,y)
plt.title(f"y = sin(2*pi*f*t) fs = {fs}")
plt.subplot(222)
plt.plot(np.linspace(-(fs/2), fs/2-1, y.shape[0]), np.fft.fftshift(abs(np.fft.fft(y))))
plt.title(f"(abs(fft(y)) fs = {fs}")

min: -0.9998507259473712 max: 1.0

Text(0.5, 1.0, '(abs(fft(y)) fs = 200')

סעיף ה: מהתוצאות ניתן לראות שכאשר דוגמים סינוס לא בכפולות שלמות, אנחנו מקבלים אות מזוגזג. כלומר, אנו דוגמים את הסינוס באפליטודה שונה בכל דגימה. כדי לקבל סינוס חלק, נדגום את האות בכפולה שלומה ובכך נדגום בכל מחזור את אותה אפליטודה של האות.

3.4

סעיף ב: הערכים שמפריעים לנו הם ערכים הגבוהים מ255, אם נוריד ערך קבוע מהמטריצה (שווה ערך לערך DC קבוע כלשהו) נוכל לטפל בבעיה ובכך למנוע ממצב בו מגיעים לערכי רוויה.

סעיף ג: הנירמול יגדיל את הטווח הדינמי ובכך לא יתרום מאחר והערכים נמצאים כבר ברוויה, במקרה הזה הטווח הדינמי מנוצל מאחר וקיימים ערכים מ0 ועד 255.

3.6

b. אורך בנוסחה נקבע לפי גודל האות שאנחנו רוצים לשכפל ועוד במרווחים כלומר אצלנו האות הוא בגודל 5 ועוד מרווח של 3 אפסים נקבל 8

f = np.array([1, 2, 3, 4, 5])
print(f"f = {f}")
size = f.shape[0]
zerosPadding = 3
size = size + zerosPadding
h = np.array([1]+[0]*size)
numberDuplicate = 4 
h = np.matlib.repmat(h,1,5).flatten()
m1 = len(f)
m2 = h.shape[0]
f = padded_zeros(f, m2+m1-1)
g = np.convolve(h,f)
print(f"h = {h}")
print("Result conv(f,h):")
print(g[:m2])

f = [1 2 3 4 5]
h = [1 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 1
 0 0 0 0 0 0 0 0]
Result conv(f,h):
[1. 2. 3. 4. 5. 0. 0. 0. 0. 1. 2. 3. 4. 5. 0. 0. 0. 0. 1. 2. 3. 4. 5. 0.
 0. 0. 0. 1. 2. 3. 4. 5. 0. 0. 0. 0. 1. 2. 3. 4. 5. 0. 0. 0. 0.]

3.5 Campbel-Robson

#%%
## create a vector of exponential values for frequencies
min_freq = 0.1 # cycles per image
N=400
min_amplitude = 1
max_amplitude = 255
max_freq = 50 # cycles per image
f = np.exp(np.linspace(np.log(min_freq), np.log(max_freq), N).reshape((1,-1)))
a = np.exp(np.linspace(np.log(min_amplitude), np.log(max_amplitude), N).reshape((-1,1)))
img = 128 + a * np.sin(2 * np.pi * f)
plt.figure(figsize=(15,15))
plt.subplot(131)
plt.title("Campbell-Robson")
plt.imshow(img, cmap='gray', vmin=0, vmax=255)
plt.subplot(132)
plt.title("Magnitude")
amp_img = np.uint8(np.matlib.repmat(a,1,img.shape[1]))
plt.imshow(amp_img, cmap='gray', vmin=0, vmax=255)
plt.subplot(133)
cos_img = uint8c(np.matlib.repmat( 255*np.cos(2 * np.pi * f)+128,img.shape[1],1))
plt.title("Cosine")
plt.imshow(cos_img, cmap='gray', vmin=0, vmax=255)
plt.figure(figsize=(10,10))
plt.subplot(211)
plt.title("Amplitude as a function of row nember")
plt.xlabel("Row")
plt.ylabel("Amp")
plt.plot(a)
plt.subplot(212)
plt.title("Frequency as a function of row nember")
plt.xlabel("Row")
plt.ylabel("Frequency")
plt.plot(f.T)

[<matplotlib.lines.Line2D at 0x7f7a4fed7590>]

fmax = 50hz fmin = 1hz AmpMax = 255 AmpMin = 1

3.7

a. LPF is used to remove high frequency noise and smooth pictures HPF is used for edge detection Median filter is most commonly used to remove salt and pepper noise b. filters are used in cascade when we want diffrent effects on the picture, for example, if we have and picture blurred filter with salt and pepper noise we would like to use median and HPF\sharpening filters in cascade c. for 2 cascade filters(MxM) on a picture of NxN the complexity would be O(M^4*N^4), we can use the properties of convolution and compute the convolution of the smaller matrices first. if N < M than we could compute the filters together and than use the product on the picture. d. a median filter cannot be used for the solution becuase it is not a filter that applied by convolution.

4 ביצוע הניסויים

4.1 - 4.5 בוצע במהלך הדוח

4.6

5 דו"ח מסכם לניסוי

5.1 - 5.5 בוצע במהלך הדוח

5.6

imfilter

הפונקציה אשר ניתן להשתמש בה באופן רב ממדי לעומת conv2 אשר רצה על דו מימד

summery

1)From the warmout excercises we concluded that a transition between two distant colors is translated to a high frequency it is easy for our eye to notice. 2) Images with diffrent noises requires specific filters, if the wrong filter is used it can make the image worse. 3) when using filters in cascade it is important to notice in what order we use the filters, for example, for a blurred image with salt and pepper noise we must use median filter first and sharpening filter after, if we switch the order the image will be substantially worse. 4) using one byte per pixel can sometimes cause problems, for example when adding images that contains sines over dc of 128 we reach values bigger than 255, converting to double and normalizing back to 0-255 can solve the problem.